function [min_scale,max_scale,mina,maxa,min_scalem,min_privatesharem, min_publicsharem,max_scalem,max_privatesharem, max_publicsharem,min_scalef,min_privatesharef, min_publicsharef,max_scalef,max_privatesharef, max_publicsharef] = get_scale(leisurem,leisuref, hworkm, hworkf, consumption,wagem,wagef,nonlabor,divorce_costsM,divorce_costsF,divorce_costsMF,ub)
%%
% LP to get bounds on a) economies of scale and the corresponding entries of A matrix, b) Male RICEB and its components and c) Female RICEB and its components
% Program uses YALMIP (https://yalmip.github.io) connected to CPLEX or GUROBI

% INPUT VARIABLES:
% leisurem[households]      -- leisure per week of the male
% leisuref[households]      -- leisure per week of the female
% hworkm[households]        -- time spend on home production by the male
% hworkf[households]        -- time spend on home production by the female
% consumption[households]   -- weekly Hicksian consumption
% wagem[households]         -- wage rate of the male
% wagef[households]         -- wage rate of the female
% nonlabor[households]      -- weekly nonlabor income of the household
% ub                        -- upper bound on economies of scale of leisure (aleisure)
% DivorceCostsM             -- Divorce Costs for individual rationality of males
% DivorceCostsF             -- Divorce Costs for individual rationality of females
% DivorceCostsMF            -- Divorce Costs for no blocking pairs 

% OUTPUT VARIABLES:
% min_scale                 -- minimum economies of scale
% max_scale                 -- maximum economies of scale
% mina                      -- 5 X 1 vector containing degree of publicness (leisure male, leisure female, domestic work male, domestic work female, Hicksian good) when economies of scale is minimized
% mina                      -- 5 X 1 vector containing degree of publicness (leisure male, leisure female, domestic work male, domestic work female, Hicksian good) when economies of scale is maximized
% min_scalem                -- minimum RICEB male
% min_privatesharem         -- share of private cosnumption when male RICEB is minimized
% min_publicsharem          -- share of public cosnumption when male RICEB is minimized
% max_scalem                -- maximum RICEB male
% max_privatesharem         -- share of private cosnumption when male RICEB is maximized
% max_publicsharem          -- share of public cosnumption when male RICEB is maximized
% min_scalef                -- minimum RICEB female
% min_privatesharef         -- share of private cosnumption when female RICEB is minimized
% min_publicsharef          -- share of public cosnumption when female RICEB is minimized
% max_scalef                -- maximum RICEB female
% max_privatesharef         -- share of private cosnumption when female RICEB is maximized
% max_publicsharef          -- share of public cosnumption when female RICEB is maximized

%%

% ==================
% defining constants
% ==================

workh = 112;                                          % total potential working hours
num_of_hh = length(leisurem);                         % total number of households in the data.


% ===========================
% defining decision variables
% ===========================

privateqM = sdpvar(num_of_hh,1);                      % private consumption of Hicksian good by male
privateqF = sdpvar(num_of_hh,1);                      % private consumption of Hicksian good by female  
privatemhworkm = sdpvar(num_of_hh,1);                 % private consumption of male's domestic work by male
privatemhworkf = sdpvar(num_of_hh,1);                 % private consumption of male's domestic work by female
privatefhworkm = sdpvar(num_of_hh,1);                 % private consumption of female's domestic work by male
privatefhworkf = sdpvar(num_of_hh,1);                 % private consumption of female's domestic work by female  
privateleisurem = sdpvar(num_of_hh,1);                % private consumption of leisure by male
privateleisuref = sdpvar(num_of_hh,1);                % private consumption of leisure by female

aleisurem = sdpvar(1);                                % degree of publicness in leisure of male
aleisuref = sdpvar(1);                                % degree of publicness in leisure of female 
ahworkm = sdpvar(1);                                  % degree of publicness in domestic work of male  
ahworkf = sdpvar(1);                                  % degree of publicness in domestic work of female 
aq = sdpvar(1);                                       % degree of publicness in Hicksian good 

nl_incomeM = sdpvar(num_of_hh,1);                     % non-labor income of male  
nl_incomeF = sdpvar(num_of_hh,1);                     % non-labor income of female

scale = sdpvar(num_of_hh,1);                          % economies of scale
scalem = sdpvar(num_of_hh,1);                         % male RICEB 
scalef = sdpvar(num_of_hh,1);                         % female RICEB
privatesharem = sdpvar(num_of_hh,1);                  % private share in male RICEB
privatesharef = sdpvar(num_of_hh,1);                  % private share in female RICEB
publicshare = sdpvar(num_of_hh,1);                    % public share in the household



%% 

% ==========================
% setting up the constraints
% ==========================

Constraints = [];

% private consumption should be balanced:
for i=1:num_of_hh
    Constraints=[Constraints, ...
        consumption(i)*(1-aq) == privateqM(i) + privateqF(i), privateqM(i) >=0, privateqF(i) >= 0,...
        hworkm(i)*(1-ahworkm) == privatemhworkm(i) + privatemhworkf(i), privatemhworkm(i) >= 0, privatemhworkf(i) >= 0,...
        hworkf(i)*(1-ahworkf) == privatefhworkm(i) + privatefhworkf(i), privatefhworkm(i) >=0 , privatefhworkf(i) >= 0,...
        leisurem(i)*(1-aleisurem) == privateleisurem(i),...
        leisuref(i)*(1-aleisuref) == privateleisuref(i)];
end


% nonlabor income should be balanced:
for i=1:num_of_hh
         Constraints = [Constraints, nonlabor(i) == nl_incomeM(i)+ nl_incomeF(i)];
end

% bounding the possible sharing of nonlabor income between partners between [.4,.6]
for i=1:num_of_hh
        if nonlabor(i)>=0
            Constraints = [Constraints,
                .4*nonlabor(i) <= nl_incomeM(i) <= .6*nonlabor(i), .4*nonlabor(i) <= nl_incomeF(i) <= .6*nonlabor(i)];
        else 
            Constraints = [Constraints,
                .6*nonlabor(i) <= nl_incomeM(i) <= .4*nonlabor(i), .6*nonlabor(i) <= nl_incomeF(i) <= .4*nonlabor(i)];
        end
end


% Individual rationality conditions:
for i=1:num_of_hh
       % IR-Male
       Constraints = [Constraints,
           (wagem(i)*workh)*divorce_costsM(i) + nl_incomeM(i) <= privateleisurem(i)*wagem(i) + privatemhworkm(i)*wagem(i) + privatefhworkm(i)*wagef(i) + privateqM(i) + hworkm(i)*ahworkm*wagem(i) + hworkf(i)*ahworkf*wagef(i) + aq*consumption(i) + leisurem(i)*aleisurem*wagem(i) + leisuref(i)*aleisuref*wagef(i)];
       % IR-Female
       Constraints = [Constraints,
           (wagef(i)*workh)*divorce_costsF(i) + nl_incomeF(i) <= privateleisuref(i)*wagef(i) + privatemhworkf(i)*wagem(i) + privatefhworkf(i)*wagef(i) + privateqF(i) + hworkm(i)*ahworkm*wagem(i) + hworkf(i)*ahworkf*wagef(i) + aq*consumption(i) + leisurem(i)*aleisurem*wagem(i) + leisuref(i)*aleisuref*wagef(i)];
end


% No blocking pairs condition: 
for i=1:num_of_hh
    for j=1:num_of_hh
        if i ~= j
           Constraints = [Constraints, 
                (wagem(i) + wagef(j))*workh*divorce_costsMF(i,j) + nl_incomeM(i) + nl_incomeF(j) <= privateleisurem(i)*wagem(i) + privatemhworkm(i)*wagem(i) + privatefhworkm(i)*wagef(j) + privateqM(i) + privateleisuref(j)*wagef(j) + privatemhworkf(j)*wagem(i) + privatefhworkf(j)*wagef(j) + privateqF(j) + max(hworkm(i),hworkm(j))*ahworkm*wagem(i) + max(hworkf(i),hworkf(j))*ahworkf*wagef(j) + aq*max(consumption(i),consumption(j)) + max(leisurem(i),leisurem(j))*aleisurem*wagem(i) + max(leisuref(i),leisuref(j))*aleisuref*wagef(j)];
        end
    end
end
   
% Constraints on degree of publicness 
Constraints = [Constraints,  0 <= aq <= 1, 0 <= aleisurem <= ub, 0 <= aleisuref <= ub, 0 <= ahworkm <= 1, 0 <= ahworkf <= 1];

% Defein economies of scale and RICEB measures
for i=1:num_of_hh
    Constraints = [Constraints,...
        scale(i) == ((1+aleisurem)*wagem(i)*leisurem(i) + (1+aleisuref)*wagef(i)*leisuref(i) + (1+ahworkm)*wagem(i)*hworkm(i)+ (1+ahworkf)*wagef(i)*hworkf(i) + (1+aq)*consumption(i))/(wagem(i)*workh + wagef(i)*workh + nonlabor(i)),...
        scalef(i) == (privateleisuref(i)*wagef(i) + privatemhworkf(i)*wagem(i) + privatefhworkf(i)*wagef(i) + privateqF(i) + aleisurem*wagem(i)*leisurem(i) +  aleisuref*wagef(i)*leisuref(i) + ahworkm*wagem(i)*hworkm(i)+ ahworkf*wagef(i)*hworkf(i) + aq*consumption(i))/(wagem(i)*workh + wagef(i)*workh + nonlabor(i)),...
        scalem(i) == (privateleisurem(i)*wagem(i) + privatemhworkm(i)*wagem(i) + privatefhworkm(i)*wagef(i) + privateqM(i) + aleisurem*wagem(i)*leisurem(i) +  aleisuref*wagef(i)*leisuref(i) + ahworkm*wagem(i)*hworkm(i)+ ahworkf*wagef(i)*hworkf(i) + aq*consumption(i))/(wagem(i)*workh + wagef(i)*workh + nonlabor(i)),...
        privatesharef(i) == (privateleisuref(i)*wagef(i) + privatemhworkf(i)*wagem(i) + privatefhworkf(i)*wagef(i) + privateqF(i))/(wagem(i)*workh + wagef(i)*workh + nonlabor(i)),...
        privatesharem(i) == (privateleisurem(i)*wagem(i) + privatemhworkm(i)*wagem(i) + privatefhworkm(i)*wagef(i) + privateqM(i))/(wagem(i)*workh + wagef(i)*workh + nonlabor(i)),...
        publicshare(i) == (aleisurem*wagem(i)*leisurem(i) +  aleisuref*wagef(i)*leisuref(i) + ahworkm*wagem(i)*hworkm(i)+ ahworkf*wagef(i)*hworkf(i) + aq*consumption(i))/(wagem(i)*workh + wagef(i)*workh + nonlabor(i))];
end


%%
% ===================
% Solving the problem
% ===================

% Specifying Objective Function
u = sdpvar(num_of_hh,1);              % vector of parameters necessasry to specify the function
Objective1 = sum(u.*scale);           % objective function that uses parameter -- u would be using as +/- unit vector 
Objective2 = sum(u.*scalem);          % objective function that uses parameter -- u would be using as +/- unit vector 
Objective3 = sum(u.*scalef);          % objective function that uses parameter -- u would be using as +/- unit vector 


%% Solving the problem
options = sdpsettings('verbose',0,'solver','gurobi');
Opt1 = optimizer(Constraints,Objective1,options,u,[scale;aleisurem;aleisuref;ahworkm;ahworkf;aq]); 
Opt2 = optimizer(Constraints,Objective2,options,u,[scalem;privatesharem;publicshare]); 
Opt3 = optimizer(Constraints,Objective3,options,u,[scalef;privatesharef;publicshare]); 

% eye matrix (diagonal 1 matrix) which would be used to solve the parametrizations of the problem
param = eye(num_of_hh,num_of_hh);

% Economies of scale
for i=1:num_of_hh
    [xvalue, errorcode] = Opt1(param(:,i));
    if errorcode~=0
        display('Something went wrong (min)');
    end
    min_scale(i) = xvalue(i);
    mina(i,1) = xvalue(num_of_hh+1);
    mina(i,2) = xvalue(num_of_hh+2);
    mina(i,3) = xvalue(num_of_hh+3);
    mina(i,4) = xvalue(num_of_hh+4);
    mina(i,5) = xvalue(num_of_hh+5);
    
    [xvalue, errorcode] = Opt1(-param(:,i));
    if errorcode~=0
        display('Something went wrong (max)');
    end
    max_scale(i) = xvalue(i);
    maxa(i,1) = xvalue(num_of_hh+1);
    maxa(i,2) = xvalue(num_of_hh+2);
    maxa(i,3) = xvalue(num_of_hh+3);
    maxa(i,4) = xvalue(num_of_hh+4);
    maxa(i,5) = xvalue(num_of_hh+5);
  end

% Male RICEB
for i=1:num_of_hh
    [xvalue, errorcode] = Opt2(param(:,i));
    if errorcode~=0
        display('Something went wrong (min)');
    end
    
    min_scalem(i) = xvalue(i);
    min_privatesharem(i) = xvalue(num_of_hh+i);
    min_publicsharem(i) = xvalue(2*num_of_hh+i);
    
    [xvalue, errorcode] = Opt2(-param(:,i));
    if errorcode~=0
        display('Something went wrong (max)');
    end
    max_scalem(i) = xvalue(i);
    max_privatesharem(i) = xvalue(num_of_hh+i);
    max_publicsharem(i) = xvalue(2*num_of_hh+i);

end

% Female RICEB
for i=1:num_of_hh
    [xvalue, errorcode] = Opt3(param(:,i));
    if errorcode~=0
        display('Something went wrong (min)');
    end
    
    min_scalef(i) = xvalue(i);
    min_privatesharef(i) = xvalue(num_of_hh+i);
    min_publicsharef(i) = xvalue(2*num_of_hh+i);

    [xvalue, errorcode] = Opt3(-param(:,i));
    if errorcode~=0
        display('Something went wrong (max)');
    end
    max_scalef(i) = xvalue(i);
    max_privatesharef(i) = xvalue(num_of_hh+i);
    max_publicsharef(i) = xvalue(2*num_of_hh+i);

end

return